open(BUFF,"whole_gprobs");
#open(BUFF,"sample_gprobs");
$line = <BUFF>;
open(WRBUFF,">final_r2_values");
open(BUFF1,"whole_r2");
while($liny = <BUFF1>)
{
	chomp($liny);
	@arry = split(/\t/,$liny);
	if($arry[1] eq "NaN")
	{
		$arry[1] = 0;
	}
	$hash{$arry[0]} = $arry[1];
}
#print "major hetero minor z w u f P0 P1 P2 dosage\n";
while($line = <BUFF>)
{
	chomp($line);
	@array = split(" ",$line);
	$rsid = $array[0];
	$f = 0;
	@z1 = ();
	@w1 = ();
	@u1 = ();
	$sumz = 0;
        $sumw = 0;
        $sumu = 0;
        $sum_zu = 0;
        $z_sqr = 0;
	
	for($i=3;$i<@array;$i++)
	{
		$major = $array[$i];
		$i++;	
		$hetero = $array[$i];
		$i++;
		$minor = $array[$i];
		$dosage = 0*$major+$hetero+$minor*2;
		
		if($dosage > 1 )
		{
			$P2=$dosage-1;
			$P1=2-$dosage;
			$P0 =0;
		}
		elsif($dosage < 1 )
		{
			$P2=0;
            		$P1 = $dosage;
           		$P0=1-$dosage;

		}
		elsif($dosage == 1)
		{
			$P0 = 0;
			$P1 = 1;
			$P2 = 0;
		}
		elsif($dosage == 0)
                {
                        $P0 = 1;
                        $P1 = 0;
                        $P2 = 0;
                }
		elsif($dosage == 2)
                {
                        $P0 = 0;
                        $P1 = 0;
                        $P2 = 2;
                }
		else
		{
			print "something wrong dosage $dosage\n";
		}
		#calculating the z value	
		if($dosage >= 1.5)
		{
			$z = 2;
		}
		elsif($dosage <= 0.5)
		{
			$z =0;
		}
		else
		{
			$z =1
		}
		#calculating u and w values
		$w = $P1+$P2*4;
		$u = $P1+$P2*2;
		$f++;
		#number
		$sumz = $sumz + $z;
		$sumw = $sumw + $w;
		$sumu = $sumu + $u;
		$sum_zu = $sum_zu + ($z*$u);
		$z_sqr = $z_sqr + ($z*$z);
		#print "test $major $hetero $minor $z $w $u $f $P0 $P1 $P2 $dosage\n";
	}
		#print "$sumz\t$sumw\t$sumu\t$sum_zu\t$z_sqr\n";
		$num = (($sum_zu -(1/$f)*($sumu*$sumz))*($sum_zu -(1/$f)*($sumu*$sumz)));
		$deno = ($sumw-(1/$f)*($sumu)*($sumu))*($z_sqr-(1/$f)*($sumz*$sumz));
		if($deno == 0)
		{
			$r2 = 0;
		}
		else
		{
			$r2 = $num/$deno;
		}
		if(exists($hash{$rsid}))
		{
			print WRBUFF $rsid."\t".$r2."\t".$hash{$rsid}."\n";
		}
		else
		{
			print "no r2 value from beagle $rsid\n";
		}
}
